        subroutine MP3(Num_ele,NORB,E_orb,U,E_mp3)
        implicit none
        integer :: i,j,k,l,a,b,c,d   !!occ:i,j,k,l vir:a,b,c,d
        integer :: NORB,occ,Num_ele
        real*8  :: U(NORB,NORB,NORB,NORB)
        real*8  :: E_orb(NORB)
        real*8  :: tem1,tem2,tem3,tem1now,E_mp3
        tem1=0.0d0;tem2=0.0d0;tem3=0.0d0;tem1now=0.0d0
        occ=Num_ele/2
        !!!part 1
        do i=1,occ
           do a=occ+1,NORB
              do j=1,occ
                 do b=occ+1,NORB
                    do c=occ+1,NORB
                       do d=occ+1,NORB
                       tem1now=4.0d0*(U(i,a,j,b)*U(a,c,b,d)*U(i,c,j,d)&
                                     +U(i,a,j,b)*U(a,d,b,c)*U(i,d,j,c)&
                                     +U(i,b,j,a)*U(a,d,b,c)*U(i,c,j,d)&
                                     +U(i,b,j,a)*U(a,c,b,d)*U(i,d,j,c))
                       tem1now=tem1now-2.0d0*&
                               (U(i,b,j,a)*U(a,c,b,d)*U(i,c,j,d)&
                               +U(i,a,j,b)*U(a,d,b,c)*U(i,c,j,d)&
                               +U(i,a,j,b)*U(a,c,b,d)*U(i,d,j,c)&
                               +U(i,b,j,a)*U(a,d,b,c)*U(i,d,j,c))
                        tem1now=tem1now/(E_orb(i)+E_orb(j)-E_orb(a)-E_orb(b))
                        tem1now=tem1now/(E_orb(i)+E_orb(j)-E_orb(c)-E_orb(d))
                        tem1=tem1+tem1now
                       enddo
                    enddo
                 enddo
              enddo
           enddo
        enddo
        !write(*,*)'tem1',tem1
        !!!part 2
        tem1now=0.0d0
        do i=1,occ
           do a=occ+1,NORB
              do j=1,occ
                 do b=occ+1,NORB
                    do k=1,occ
                       do l=1,occ
                       tem1now=4.0d0*(U(i,a,j,b)*U(i,k,j,l)*U(k,a,l,b)&
                                     +U(i,a,j,b)*U(i,l,j,k)*U(k,b,l,a)&
                                     +U(i,b,j,a)*U(i,l,k,j)*U(k,a,l,b)&
                                     +U(i,b,j,a)*U(i,k,j,l)*U(k,b,l,a))
                       tem1now=tem1now-2.0d0*&
                                (U(i,a,j,b)*U(i,l,j,k)*U(k,a,l,b)&
                                +U(i,a,j,b)*U(i,k,j,l)*U(k,b,l,a)&
                                +U(i,b,j,a)*U(i,k,j,l)*U(k,a,l,b)&
                                +U(i,b,j,a)*U(i,l,j,k)*U(k,b,l,a)) 
                       tem1now=tem1now/(E_orb(i)+E_orb(j)-E_orb(a)-E_orb(b))
                       tem1now=tem1now/(E_orb(k)+E_orb(l)-E_orb(a)-E_orb(b))
                       tem2=tem2+tem1now
                       enddo
                    enddo
                 enddo
              enddo
           enddo
        enddo
        !write(*,*)'tem2',tem2
        tem1now=0.0d0 
        do i=1,occ
           do a=occ+1,NORB
              do j=1,occ
                 do b=occ+1,NORB
                    do c=occ+1,NORB
                       do k=1,occ
                        tem1now=8d0*(U(i,a,j,b)*U(j,b,c,k)*U(i,a,k,c))
                        tem1now=tem1now-4d0*&
                                    (U(i,b,j,a)*U(j,b,c,k)*U(i,a,k,c)&
                                    +U(i,a,j,b)*U(j,k,c,b)*U(i,a,k,c)&
                                    +U(i,a,j,b)*U(j,b,c,k)*U(i,c,k,a)&
                                    +U(i,b,j,a)*U(j,k,c,b)*U(i,c,k,a))
                        tem1now=tem1now+2d0*&
                                (U(i,b,j,a)*U(j,k,c,b)*U(i,a,k,c)&
                                +U(i,b,j,a)*U(j,b,c,k)*U(i,c,k,a)&
                                +U(i,a,j,b)*U(j,k,c,b)*U(i,c,k,a))
                        tem1now=tem1now/(E_orb(i)+E_orb(j)-E_orb(a)-E_orb(b))
                        tem1now=tem1now/(E_orb(i)+E_orb(k)-E_orb(c)-E_orb(a))
                        tem3=tem3+tem1now
                       end do
                    end do
                 end do
              end do
            end do
        end do

        E_mp3=0.125*tem1+0.125*tem2+tem3
        endsubroutine MP3       
         
